QK Box 8.2

Loyn (1987) selected 56 forest patches in southeastern Victoria, Australia, and related the abundance of forest birds in each patch to six predictor variables: patch area (ha), distance to nearest patch (km), distance to nearest larger patch (km), grazing stock (1 to 5 indicating light to heavy), altitude (m) and years since isolation (years). The aim was to develop a predictive model relating bird abundance to these predictors and to assess which predictors were most important for bird abundance.

Tawny Frogmouth, Podargus strigoides. Mick Keough, CC BY 4.0

Eastern Yellow Robin, Eopsaltria australis. Emmet Keough, CC BY 4.0

This data set was used in the first edition and is available here

Preliminaries

First, load the required packages (car, lm.beta, hier.part)

Import loyn data file (loyn.csv)

loyn <- read.csv("../data/loyn.csv")
head(loyn,10)

NA

Check distributions of predictors and look for correlations

We will use the scatterplotMatrix function from the car package

scatterplotMatrix(~abund+area+dist+graze+alt+yearisol, data=loyn, cex=0.33, regLine=FALSE, diagonal=list(method='boxplot'))

abund looks symmetrical but two unusual observations for area resulting in non-linear relationships with abund and one outlier for distance to nearest patch

Try log(10) transformation of area and dist

scatterplotMatrix(~abund+log10(area)+log10(dist)+graze+alt+yearisol, data=loyn, cex=0.33, regLine=FALSE, diagonal=list(method='boxplot'))

Plots look much better. Now transform the variables in loyn

loyn$logarea <- log10(loyn$area)
loyn$logdist <- log10(loyn$dist)

Now check for collinearity among predictors (using larea and ldist)

Use correlations and VIF

scatterplotMatrix(~logarea+logdist+graze+alt+yearisol, data=loyn, cex=0.33, regLine=FALSE, diagonal=list(method='boxplot'))

cor(loyn[,c('logarea','logdist','graze','alt','yearisol')])

            logarea     logdist      graze        alt    yearisol
logarea   1.0000000  0.30216662 -0.5590886  0.2751428  0.27841452
logdist   0.3021666  1.00000000 -0.1426392 -0.2190070 -0.01957223
graze    -0.5590886 -0.14263922  1.0000000 -0.4071671 -0.63556710
alt       0.2751428 -0.21900701 -0.4071671  1.0000000  0.23271541
yearisol  0.2784145 -0.01957223 -0.6355671  0.2327154  1.00000000

vif(lm(abund~ logarea+logdist+graze+alt+yearisol, data=loyn))

 logarea  logdist    graze      alt yearisol 
1.619373 1.265608 2.522813 1.364838 1.735318 

Collinearity OK so now fit linear model relating response to predictor

loyn.lm <- lm(abund~ logarea+logdist+graze+alt+yearisol, data=loyn)

Examine regression diagnostics (residual plot and Cooks D)

plot(loyn.lm)

augment(loyn.lm)

Display results of model fitting

Get parameter estimates, their confidence intervals, and tests

tidy(loyn.lm, conf.int=TRUE)

Get standarised regression coefficients

Use lm.beta function from lm.beta package

lm.beta.loyn <- lm.beta(loyn.lm)
tidy(lm.beta.loyn, conf.int=TRUE)

Show the partial regression (added-variable) plots

avPlots(loyn.lm, ask=F)

Run hierarchical partitioning (used later)

We need two dataframes that are subsets of loyn containing the response variable and the predictors

loyn_abund<-loyn$abund
loyn_pred<-subset(loyn, select = c("logarea","logdist","alt","yearisol","graze"))

hier.part(loyn_abund, loyn_pred, family="gaussian", gof="Rsqu")

$gfs
 [1] 0.0000000 0.5476530 0.0160588 0.1488696 0.2533690 0.4658218 0.5579841 0.5835769 0.6434808 0.6527344 0.1957330 0.2720289 0.4667023 0.3297010 0.4797883 0.4739432 0.5852970
[18] 0.6479714 0.6609983 0.6628587 0.6610593 0.6732151 0.3707161 0.4845809 0.4758029 0.4887284 0.6634401 0.6651442 0.6787941 0.6823418 0.4960985 0.6843360

$IJ

$I.perc

$params
$params$full.model
[1] "y ~ logarea + logdist + alt + yearisol + graze"

$params$family
[1] "gaussian"

$params$link
[1] "default"

$params$gof
[1] "Rsqu"